home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Precision Software Appli…tions Silver Collection 4
/
Precision Software Applications Silver Collection Volume 4 (1993).iso
/
stats
/
yamp15.exe
/
VIRTDISK.CPP
< prev
next >
Wrap
C/C++ Source or Header
|
1992-08-22
|
28KB
|
1,197 lines
/*
YAMP - Yet Another Matrix Program
Version: 1.5
Author: Mark Von Tress, Ph.D.
Date: 08/28/92
Copyright(c) Mark Von Tress 1992
Portions of this code are (c) 1991 by Allen I. Holub and are used by
permission
DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
FROM USE OF THIS PROGRAM.
*/
// #include "virt.h"
#include "alloc.h"
//////////////////////////////////////////// string functions
strtype::strtype(strtype &str) // copy constructor
{ int len = MAXCHARS;
len = strlen(str.name);
len = (len >= MAXCHARS) ? MAXCHARS - 2 : len;
strncpy(name, str.name, len);
name[len] = '\0';
}
strtype::strtype(char *str)
{ int len = MAXCHARS;
len = strlen(str);
len = (len >= MAXCHARS) ? MAXCHARS - 2 : len;
strncpy(name, str, len);
name[len] = '\0';
}
strtype::strtype(void)
{
name[0] = '\0';
}
strtype strtype::operator + (strtype str)
{
int len1, len2;
len1 = strlen(name);
len2 = strlen(str.name);
int len =(len2 + len1 >= MAXCHARS) ? MAXCHARS - len1 - 3 : len2;
len = (len < 0 || len >= MAXCHARS - 1) ? 0 : len;
strncat(name, str.name, len);
name[len + len1] = '\0';
return *this;
}
strtype strtype::operator + (const char *str)
{
int len1 = strlen(name);
int len2 = strlen(str);
int len = (len2 + len1 >= MAXCHARS) ? MAXCHARS - len1 - 3 : len2;
len = (len < 0 || len >= MAXCHARS - 1) ? 0 : len;
len = (len >= MAXCHARS) ? 0 : len;
strncat(name, str, len);
name[len + len1] = '\0';
return *this;
}
strtype strtype::operator = (strtype str)
{
int len = (strlen(str.name) >= MAXCHARS) ? MAXCHARS - 3
: strlen(str.name);
strncpy(name, str.name, len);
name[len] = '\0';
return *this;
}
strtype strtype::operator = (const char *str)
{
int len = (strlen(str) >= MAXCHARS) ? MAXCHARS - 2 : strlen(str);
strncpy(name, str, len);
name[len] = '\0';
return *this;
}
//////////////////////////////////// virtual memory allocation functions
// Virtual memory routines used by permission
// Portions of this code are (c) 1991 by Allen I. Holub and are used by
// permission
// The code is found in Programmer's Journal volume 8.5. (1990)
//
//#define DEBUG 1
#ifdef DEBUG
#define D(x) x
#define ND(x)
#else
#define D(x)
#define ND(x) x
#endif
#define READ 1
#define WRITE 0
#define TNAME "$$$$vmem.tmp"
#define BLK_SIZE 512
#define SIGNATURE 0x5a5a
static char Vname[128];
static unsigned Filesize = 0;
static int Vfd = -1;
static int nobj = 0;
static hdr *Freelist = NULL;
static hdr *new_hdr(unsigned nblocks, unsigned offset);
static unsigned enlarge(unsigned blocks_reqd);
static hdr *add_vmem(long num_ele, int ele_size, hdr *header);
static int load_block(hdr *p, unsigned req_block);
static int block_io(int doread, unsigned block, char *buf);
/*--------------------------------------------*/
void pheader(hdr *p)
{
char junk[12] = " <-Freelist";
char s[13] = " ";
if (!p) printf("NULL\n");
else {
if (p == Freelist) strcpy(s, junk);
printf("hdr 0x%lX :active %1d offset = %u, nblocks = %u, next=0x%lX %s\n",
p, p->active, p-> offset, p-> nblocks, p-> next.h, s);
}
}
void pfreelist(void)
{
hdr *p = NULL;
if (!(p = Freelist)) printf("Freelist is empty.\n");
else {
printf("Freelist is:\n");
do {
pheader(p);
} while ((p = p->next.h) != Freelist);
}
}
/*-------------------------------------*/
int vopen(void)
{
char *env;
int len;
if (!(env = getenv("TMP")))
strcpy(Vname, TNAME);
else {
if ((len = strlen(env)) > sizeof (Vname) - 16) {
errno = E2BIG;
return 0;
}
sprintf(Vname,
(len >= 2 && env[len - 2] != '/' && env[len - 2] != '\\')?
"%s/%s" : "%s%s", env, TNAME);
}
if ((Vfd = creat(Vname, S_IREAD | S_IWRITE)) == -1) return 0;
close(Vfd);
if ((Vfd = open(Vname, O_RDWR | O_BINARY)) == -1) return 0;
return 1;
}
/*---------------------------------*/
int vclose(void)
{
hdr *p = NULL, *newp = NULL;
if (Vfd == -1) {
errno = EBADF;
return 0;
}
else {
if (close(Vfd) == -1) return 0;
Vfd = -1;
ND(if (unlink(Vname) == -1) return 0;)
if (Freelist) {
p = Freelist;
do {
newp = p->next.h;
free(p);
} while ((p = newp) != Freelist);
Freelist = NULL;
}
}
return 1;
}
/*----------------------------------*/
hdr *vmalloc(long num_ele, size_t ele_size)
{
hdr *prev, *cur, *newhdr, *eof_block;
unsigned ele_per_block;
unsigned blocks_reqd;
D(printf("Entering vmalloc(%ld,%u). ", num_ele, ele_size);)
D(pfreelist();)
ele_per_block = BLK_SIZE / ele_size;
blocks_reqd = num_ele / ele_per_block + (num_ele % ele_per_block != 0);
if (Freelist) {
prev = Freelist;
cur = Freelist->next.h;
eof_block = NULL;
while (1) {
do {
if (cur-> offset + cur-> nblocks == Filesize) eof_block = cur;
if (cur-> nblocks >= blocks_reqd) {
Freelist = prev;
if (cur-> nblocks > blocks_reqd) {
cur->nblocks -= blocks_reqd;
if (newhdr = new_hdr(blocks_reqd, cur->offset + cur->nblocks))
newhdr = add_vmem(num_ele, ele_size, newhdr);
goto end;
}
else {
if (cur-> next.h == cur) Freelist = NULL;
else prev->next.h = cur->next.h;
newhdr = add_vmem(num_ele, ele_size, cur);
goto end;
}
}
prev = cur;
cur = cur->next.h;
} while (prev != Freelist);
if (!eof_block) break;
else if (!enlarge(blocks_reqd - eof_block->nblocks)) {
newhdr = NULL;
goto end;
}
else eof_block->nblocks = blocks_reqd;
} /* end while(1) */
} /* end if(freelist) */
if (!(newhdr = new_hdr(blocks_reqd, Filesize))) goto end;
if (!enlarge(blocks_reqd)) {
free(newhdr);
newhdr = NULL;
}
newhdr = add_vmem(num_ele, ele_size, newhdr);
end :
if (newhdr != NULL) {
newhdr->active = 1;
newhdr->nrefs = 1;
}
D(printf("leaving vmalloc(%ld,%u), Return:\n", num_ele, ele_size);)
D(pheader(newhdr);)
D(pfreelist();)
D(printf("-----------------------------\n");)
return newhdr;
}
/*---------------------------------------*/
hdr *new_hdr(unsigned nblocks, unsigned offset)
{
hdr *newhdr = NULL;
if (!(newhdr = (hdr *) malloc(sizeof (hdr)))) {
errno = ENOMEM;
return NULL;
}
newhdr->nblocks = nblocks;
newhdr->offset = offset;
return newhdr;
}
/*---------------------------------------*/
unsigned enlarge(unsigned blocks_reqd)
{
unsigned size_in_blocks;
long real_size, position;
position = ((long) blocks_reqd) * ((long) BLK_SIZE) - 1L;
if ((real_size = lseek(Vfd, position, SEEK_END)) == -1) {
errno = ENOMEM;
return 0L;
}
if (write(Vfd, "", 1) == -1)
return 0L;
++ real_size;
return ((real_size / BLK_SIZE > (unsigned) ~0) ? 0 :
(Filesize = (unsigned) (real_size / BLK_SIZE)));
}
/*---------------------------------------*/
hdr *add_vmem(long num_ele, int ele_size, hdr *header)
{
vmem *p = NULL;
if (!(p = (vmem *) malloc(sizeof (vmem)))) {
errno = ENOMEM;
return NULL;
}
else if (!(p->buf = (char *) malloc(BLK_SIZE))) {
errno = ENOMEM;
return NULL;
}
else {
header->next.v = p;
p->signature = SIGNATURE;
p->dirty = 0;
p->ele_size = ele_size;
p->ele_per_blk = BLK_SIZE / ele_size;
p->num_ele = num_ele;
p->cblock = 0;
return header;
}
}
/*-------------------------------------*/
int vfree(hdr *p)
{
int rval = 1;
hdr *cur, *h;
vmem *v;
v = p->next.v;
D(printf("Entering vfree(0x%lX).\n", p);)
D(printf("disposing: "); pheader(p);)
D(pfreelist();)
if (Vfd == -1) { /* VMS system is not open */
D(printf("VFREE: VMS system closed\n");)
goto end;
}
p->nrefs = 0;
if (v-> signature == SIGNATURE) v-> signature = 0;
else {
errno = EBADF;
rval = 0;
goto end;
}
p->active = 0;
free(v->buf);
free(v);
if (!Freelist) {
Freelist = p->next.h = p;
goto end;
}
for (cur = Freelist; 1; cur = cur-> next.h) {
if ((cur->offset < p->offset && p->offset < cur->next.h->offset)
|| (cur-> offset >= cur-> next.h->offset &&
(cur-> offset < p-> offset ||
p->offset < cur-> next.h->offset)))
break;
}
if (cur == cur-> next.h) {
if (p-> offset + p-> nblocks == cur-> offset) {
cur->offset = p-> offset;
cur->nblocks += p-> nblocks;
if (p == Freelist) Freelist = cur;
free(p);
goto end;
}
else if (cur-> offset + cur-> nblocks == p-> offset) {
cur->nblocks += p-> nblocks;
if (p == Freelist) Freelist = cur;
free(p);
goto end;
}
}
if (p-> offset + p-> nblocks == cur-> next.h->offset) {
p->nblocks += cur-> next.h->nblocks;
p->next.h = cur->next.h->next.h;
if (cur-> next.h == Freelist) Freelist = cur;
free(cur->next.h);
}
else
p->next.h = cur->next.h;
if (cur-> offset + cur-> nblocks == p-> offset) {
cur->nblocks += p-> nblocks;
cur->next.h = p->next.h;
if (p == Freelist) Freelist = cur;
free(p);
}
else cur->next.h = p;
end :
D(printf("leaving vfree( 0x%lX). Returning %d. ", p, rval);)
D(pfreelist();)
D(printf("-----------------------------\n");)
return rval;
}
/*-----------------------------------------------*/
int block_io( int doread, unsigned block_num, char *buf)
{
size_t got;
long here;
here = lseek(Vfd, ((long)block_num * (long)BLK_SIZE) , SEEK_SET);
if (here == -1 ) return 0;
got = doread ? read (Vfd, buf, (size_t) BLK_SIZE)
: write(Vfd, buf, (size_t) BLK_SIZE);
if ( got == -1 ) return 0;
return (doread ? 1 : (got == BLK_SIZE));
}/*----------------------------------------------*/
int load_block(hdr *p, unsigned req_block)
{
vmem *v;
v = p->next.v;
if (req_block != v-> cblock) {
if (v-> dirty && !block_io(WRITE, p->offset + v->cblock, v->buf))
return 0;
if (!block_io(READ, p->offset + req_block, v->buf)) return 0;
v->dirty = 0;
v->cblock = req_block;
}
return 1;
}
/*-----------------------------------------*/
void *vread(hdr *p, long index)
{
unsigned block_num;
vmem *v;
v = p->next.v;
if (index < 0 || index > v-> num_ele) {
errno = ERANGE;
return 0;
}
if (v-> signature != SIGNATURE || Vfd == -1) {
errno = EBADF;
return 0;
}
block_num = index / v->ele_per_blk;
if (!load_block(p, block_num)) return NULL;
index -= block_num * v->ele_per_blk;
index *= v->ele_size;
return (v-> buf + index);
}
/*---------------------------------------------*/
void *vwrite(hdr *p, long index, char *thiss)
{
int i;
char *bp = NULL;
void *startp = NULL;
vmem *v;
v = p->next.v;
if (startp = vread(p, index)) {
bp = (char *) startp;
for (i = v-> ele_size;-- i >= 0;) *bp++ = *thiss++;
v->dirty = 1;
}
return startp;
}
/*--------------------------------------------------*/
long vele(hdr *p) { return p->next.v->num_ele; }
void vdirty(hdr *p) { p->next.v->dirty = 1; }
/*--------------------------------------------------*/
///////////////////////////////////////////// virtual double array
static double val(vdoub &v)
{
if (v.cur_ele) return *v.cur_ele;
cerr << "VMS: no previous index on vdoub (val)\n";
return 0;
}
vdoub::vdoub(long array_size)
{
if (nobj++ == 0) {
if (!vopen()) {
perror("VMS: unable to open virtural memory file");
exit(1);
}
}
cur_ele = NULL;
cur_index = 0;
signature = SIGNATURE;
if (!(hdr = vmalloc(array_size, sizeof (double)))) {
perror("VMS allocation error 1");
exit(1);
}
}
vdoub::vdoub(void)
{
if (nobj++ == 0) {
if (!vopen()) {
perror("VMS: unable to open virtural memory file");
exit(1);
}
}
cur_ele = NULL;
cur_index = 0;
signature = SIGNATURE;
if (!(hdr = vmalloc(1L, sizeof (double)))) {
Dispatch->PrintStack();
perror("VMS allocation error 2");
exit(1);
}
}
vdoub::vdoub(vdoub &src)
{
if (nobj++ == 0) {
if (!vopen()) {
perror("VMS: unable to open virtural memory file");
exit(1);
}
}
double *p;
long numele = vele(src.hdr);
signature = SIGNATURE;
cur_ele = NULL;
cur_index = 0;
if (!(hdr = vmalloc(numele, sizeof (double)))) {
perror(" VMS copy allocation error");
exit(1);
}
while (-- numele >= 0) {
if (!(p = (double *) vread(src.hdr, numele))) {
perror("VMS copy-access( read ) error");
exit(1);
}
if (!vwrite(hdr, numele, (char *) p)) {
perror("VMS copy-access( write ) error");
exit(1);
}
}
}
vdoub::~vdoub(void)
{
if (0 ==-- (hdr->nrefs)) {
if (!vfree(hdr)) {
perror("VMS deallocation error ");
exit(1);
}
}
signature = 0; // make sure to change signature in destructor
if (-- nobj == 0) vclose();
}
double vdoub::v(long index)
{
if (!(cur_ele = (double *) vread(hdr, index))) {
cerr << "hdr, index: " << hdr << " " << index;
perror("VMS selection error 1");
}
cur_index = index;
return *cur_ele;
}
vdoub &vdoub::operator[](long index)
{
if (!(cur_ele = (double *) vread(hdr, index))) {
cerr << "hdr, index: " << hdr << " " << index << " ";
perror("VMS selection error 2");
pfreelist();
exit(1);
}
cur_index = index;
return *this;
}
double vdoub::operator = (double d)
{
if (Garbage()) {
perror("VMS assignment error, source is Garbage(=)");
exit(1);
}
*cur_ele = d;
vdirty(hdr);
return d;
}
double vdoub::operator = (vdoub &v)
// copy vector v to vector t
{
if (v.Garbage()) {
perror("VMS copy error, source is Garbage");
exit(1);
}
if (!Garbage()) vfree(hdr); // replace hdr if t is active
double *p;
long numele = vele(v.hdr);
signature = SIGNATURE;
*cur_ele = *v.cur_ele;
cur_index = v.cur_index;
if (!(hdr = vmalloc(numele, sizeof (double)))) {
perror(" VMS copy allocation error");
exit(1);
}
while (-- numele >= 0) {
if (!(p = (double *) vread(v.hdr, numele))) {
perror("VMS copy-access( read ) error");
exit(1);
}
if (!vwrite(hdr, numele, (char *) p)) {
perror("VMS copy-access( write ) error");
exit(1);
}
}
return *cur_ele;
}
//////////////////////////////////////////////////matrix stack functions
MStack::MStack(void)
{
static int cnter = 0;
next = NULL;
stackloc = 0;
level = 0;
if (!cnter) {
VMatrix::Nameit("DISPATCH");
cnter = 1;
level = 1;
}
}
void MStack::Inclevel(void)
{
Dispatch->level++;
}
void MStack::Declevel(void)
{
(Dispatch-> level)--;
if (Dispatch-> level < 0) {
printf("ERROR: LEVEL has been decremented too often\n");
exit(1);
}
}
void VMatrix::NewReference(VMatrix &ROp)
{
cur_ele = NULL;
cur_index = 0;
signature = SIGNATURE;
r = ROp.r;
c = ROp.c;
// release the current header and share it with ROp.hdr
PurgeVectors();
hdr = ROp.hdr;
hdr->nrefs += 1;
}
void MStack::Push(VMatrix &ROp)
{
ROp.Garbage("Push");
D(printf("Pushing: "); Dispatch->PrintStack();)
MStack *temp = new MStack;
temp->Nameit(ROp.Getname(ROp));
temp->NewReference(ROp);
temp->next = Dispatch-> next;
Dispatch->next = temp;
temp->level = Dispatch-> level;
temp->stackloc =++ (Dispatch->stackloc);
}
VMatrix& MStack::Pop(void)
{
D(printf("Popping: "); Dispatch->PrintStack();)
if (Dispatch-> next && Dispatch-> stackloc) {
MStack *t = Dispatch->next;
Dispatch->next = Dispatch-> next-> next;
delete t;
(Dispatch-> stackloc)--;
}
else Nrerror("POP: Stack is empty, cant pop any more\n");
return *Dispatch;
}
void MStack::Cleanstack(void)
{
while (Dispatch-> next-> level >= Dispatch-> level
&& Dispatch->next-> next)
Dispatch->Pop();
}
void MStack::PrintStack(void)
{
MStack *temp = Dispatch;
int t = 1 + Dispatch->stackloc;
while (t--) {
printf("\n Stack matrix :location level r c:%d %d %d %d: name: ",
t, temp->level, temp-> r, temp-> c);
temp->Showname();
temp = temp->next;
}
printf("\n\n");
}
//////////////////////////////////////////////////////matrix class
VMatrix::VMatrix(void)
{
r = c = 1;
curvecind = 0;
Nameit("t");
signature = SIGNATURE;
// for( int i=1; i<=r; i++ ){
// for(int j=1; j<=c; j++) M(i,j) = 0;
// }
}
VMatrix::VMatrix(const char *str, int rr, int cc) :
vdoub(((long) ((long) rr) *((long) cc)))
{
r = rr;
c = cc;
curvecind = 0;
Nameit(str);
signature = SIGNATURE;
// for( int i=1; i<=r; i++ ){
// for(int j=1; j<=c; j++) M(i,j) = 0;
// }
}
VMatrix::VMatrix(VMatrix &ROp) // copy constructor
{
ROp.Garbage("Copy Constructor");
r = ROp.r;
c = ROp.c;
name = ROp.name;
curvecind = 0;
signature = SIGNATURE;
PurgeVectors(); // note a hdr is constructed in the
// vdoub constructor so it must be deleted
// first.
SetupVectors(r, c);
for (int i = 1; i <= r;++ i) {
for (int j = 1; j <= c;++ j) M(i, j) = ROp.m(i, j);
}
Dispatch->Cleanstack();
}
VMatrix::~VMatrix(void)
{
signature = 0;
}
//////////////////////////////////// matrix member functions
double VMatrix::Nrerror(const char *errormsg)
{
double x;
printf("\nMATRIX ERROR: %s \nexiting to system\n", errormsg);
x = 2;
exit(1);
return x;
}
void VMatrix::Garbage(const char *errormsg) {
#ifndef NO_CHECKING
int errorint = 0;
if (signature != SIGNATURE) errorint++;
if (errorint) {
Showname();
printf("\nFunction name: %s", errormsg);
Nrerror("Operating on Garbage");
}
#endif
return;
}
void VMatrix::SetupVectors(int rr, int cc)
{
cur_ele = NULL;
cur_index = 0;
signature = SIGNATURE;
r = rr;
c = cc;
if (!(hdr = vmalloc(((long) ((long) r) *((long) c)), sizeof (double)))){
perror("VMS allocation error 3");
exit(1);
}
}
void VMatrix::PurgeVectors(void)
{
// frees the virtual vector
Garbage("PurgeVectors");
if (!vfree(hdr)) {
perror("VMS deallocation error ");
exit(1);
}
}
void VMatrix::DisplayMat(void)
{
Garbage("DisplayMat");
int wid = Getwid();
int dec = Getdec();
printf("\n");
name.Showname();
printf("\n\n");
for (int i = 1; i <= r;++ i) {
for (int j = 1; j <= c;++ j) {
printf("%*.*lf ", wid, dec, m(i, j));
}
printf("\n");
}
printf("\n");
}
void VMatrix::InfoMat(void)
{
Garbage("InfoMat");
printf("\n");
name.Showname();
printf("\nr c: %d %d\n", r, c);
}
void VMatrix::LoadMat(void)
{
// interactive matrix loading;
char buffer[256] = { '\0' };
Garbage("LoadMat");
printf(" Enter name: ");
cin >> buffer;
Nameit(buffer);
double d;
for (int j = 1; j <= r;++ j) {
for (int k = 1; k <= c;++ k) {
cout << "Enter (" << j << "," << k << "): ";
cin >> d;
M(j, k) = d;
}
}
}
void VMatrix::Replace(VMatrix &ROp)
{
ROp.Garbage("Replace");
if (r != ROp.r || c != ROp.c) {
PurgeVectors();
SetupVectors(ROp.r, ROp.c);
}
name = ROp.name;
for (int i = 1; i <= r;++ i) {
for (int j = 1; j <= c;++ j) { M(i, j) = ROp.m(i, j);
}
}
}
void VMatrix::operator = (VMatrix &ROp)
{
Garbage("operator =");
ROp.Garbage("operator =");
D(printf("Equals "); Showname(); Dispatch->PrintStack();)
Replace(ROp);
Dispatch->Cleanstack();
}
////////////////////// end matrix member functions
/////////////////// freind functions of matrix class
///------------- addition
VMatrix& operator + (VMatrix &LOp, VMatrix &ROp)
{
LOp.Garbage("operator +");
ROp.Garbage("operator +");
if (LOp.r == ROp.r && LOp.c == ROp.c) {
VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
for (int i = 1; i <= LOp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = LOp.m(i, j) + ROp.m(i, j);
}
}
temp->name = temp-> name + LOp.name + "+" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
}
else ROp.Nrerror("Mismatched Matrix sizes in addition function\n");
return Dispatch->ReturnMat();
}
VMatrix& operator +(double scalar, VMatrix &ROp)
{
ROp.Garbage("operator s+M");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = scalar + ROp.m(i, j);
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + string + "+" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& operator +(VMatrix &ROp, double scalar)
{
ROp.Garbage("operator M+s");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = scalar + ROp.m(i, j);
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + ROp.name + "+" + string + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
////-------------subtraction
VMatrix& operator -(VMatrix &LOp, VMatrix &ROp)
{
LOp.Garbage("operator M-M");
ROp.Garbage("operator M-M");
if (LOp.r == ROp.r && LOp.c == ROp.c) {
VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
for (int i = 1; i <= LOp.r;++ i) {
for (int j = 1; j <= LOp.c;++ j) {
temp->M(i, j) = LOp.m(i, j) - ROp.m(i, j);
}
}
temp->name = temp-> name + LOp.name + "-" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
}
else Dispatch->Nrerror("Mismatched VMatrix sizes in subtraction function\n");
return Dispatch->ReturnMat();
}
VMatrix& operator -(double scalar, VMatrix &ROp)
{
ROp.Garbage("operator s-M");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = scalar - ROp.m(i, j);
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + string + "-" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& operator -(VMatrix &ROp, double scalar)
{
ROp.Garbage("operator M-s");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = ROp.m(i, j) - scalar;
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + ROp.name + "-" + string + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
//------------- unary minus
VMatrix& operator -(VMatrix &ROp)
{
ROp.Garbage("operator -");
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = -ROp.m(i, j);
}
}
temp->name = temp-> name + "-" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
//-------------- multiplication
VMatrix& operator* (VMatrix &LOp, VMatrix &ROp)
{
long double sum = 0;
LOp.Garbage("operator M*M");
ROp.Garbage("operator M*M");
if (LOp.c == ROp.r) {
VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
for (int i = 1; i <= LOp.r; i++) {
for (int j = 1; j <= ROp.c; j++) {
sum = 0.0;
for (int k = 1; k <= LOp.c; k++)
sum +=(long double)
((long double) LOp.m(i, k)) *((long double) ROp.m(k,j));
temp->M(i, j) = (double) sum;
}
}
temp->name = temp-> name + LOp.name + "*" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
}
else Dispatch->Nrerror(
"Mismatched VMatrix sizes in multiplication function\n");
return Dispatch->ReturnMat();
}
VMatrix& operator*(double scalar, VMatrix &ROp)
{
ROp.Garbage("operator s*M");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = scalar * ROp.m(i, j);
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + string + "*" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& operator* (VMatrix &ROp, double scalar)
{
ROp.Garbage("operator M*s");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = scalar * ROp.m(i, j);
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + ROp.name + "*" + string + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
//-------- elementwise multiplication
VMatrix& operator %(VMatrix &LOp, VMatrix &ROp)
{
LOp.Garbage("operator M%M");
ROp.Garbage("operator M%M");
if (LOp.r == ROp.r && LOp.c == ROp.c) {
VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
for (int i = 1; i <= LOp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j)
temp->M(i, j) = LOp.m(i, j) *ROp.m(i, j);
}
temp->name = temp-> name + LOp.name + "%" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
}
else Dispatch->Nrerror(
"Mismatched Matrix sizes in elementwise multiplication\n");
return Dispatch->ReturnMat();
}
//------------- division
VMatrix& operator /(VMatrix &LOp, VMatrix &ROp)
{
LOp.Garbage("operator M/M");
ROp.Garbage("operator M/M");
if (LOp.r == ROp.r && LOp.c == ROp.c) {
VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
for (int i = 1; i <= LOp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
double d = ROp.m(i, j);
double L = LOp.m(i, j);
temp->M(i, j) = (fabs(d) > ZERO || fabs((d - L) / d) < 1.0E-5 )
? L / d
: ROp.Nrerror(" zero divisor in elementwise division");
}
}
temp->name = temp-> name + LOp.name + "/" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
}
else Dispatch->Nrerror(
"Mismatched Matrix sizes in elementwise division\n");
return Dispatch->ReturnMat();
}
VMatrix& operator /(VMatrix &ROp, double scalar)
{
ROp.Garbage("operator M/s");
strtype string;
char buffer[MAXCHARS] = { '\0' };
if (fabs(scalar) < ZERO)
ROp.Nrerror(" zero divisor in elementwise division");
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = ROp.m(i, j) / scalar;
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + ROp.name + "/" + string + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
VMatrix& operator /(double scalar, VMatrix &ROp)
{
ROp.Garbage("operator s/M");
strtype string;
char buffer[MAXCHARS] = { '\0' };
VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
for (int i = 1; i <= ROp.r;++ i) {
for (int j = 1; j <= ROp.c;++ j) {
temp->M(i, j) = (ROp.m(i, j) != 0.0) ? scalar / ROp.m(i, j)
: ROp.Nrerror("zero divisor in scalar divison");
}
}
gcvt(scalar, NDECS, buffer);
string = buffer;
temp->name = temp-> name + string + "/" + ROp.name + ")";
Dispatch->Push(*temp);
delete temp;
return Dispatch->ReturnMat();
}
//////////////////////end of friend functions
////////////////////// matrix functions
////////////// utility functions
void VMatrix::Writeb(char *fid, VMatrix &mat)
/* write a matrix in an binary file */
{
int i;
FILE *file;
double *x;
char name[MAXCHARS] = { '\0' };
file = fopen(fid, "wb");
if (!file) Dispatch->Nrerror("WRITEB: unable to open file");
mat.Garbage("Writeb");
strncpy(name, mat.StringAddress(), 79);
i = strlen(name);
fwrite(&i, sizeof (int), 1, file);
fputs(name, file);
i = mat.r;
fwrite(&i, sizeof (int), 1, file);
i = mat.c;
fwrite(&i, sizeof (int), 1, file);
unsigned blocks = mat.hdr-> nblocks;
char *locbuf = mat.hdr-> next.v->buf;
unsigned block_num;
for (i = 0; i < blocks; i++) {
if (!load_block(mat.hdr, (unsigned) i))
Nrerror("WRITEB: could not load a virtual block");
fwrite(locbuf, sizeof (char), BLK_SIZE, file);
}
fclose(file);
mat.M(1, 1); // reset matrix to base pointer
} /* writeb */